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Abstract. Magnetization transport in a one-dimensional isotropic spin 1/2 
Heisenberg model is studied. It is shown that in a nonequilibrium steady state at 
high temperature and constant small driving the magnetization current depends on 
the system length L as '-^ l/L"-^, meaning that the diffusion constant diverges as 
^ L*^'^. Spectral properties of a superoperator governing the relaxation towards a 
nonequilibrium steady state are also discussed. 
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1. Introduction 

The Heisenberg model of nearest-neiglibor coupled spins is of high interest in theoretical 
as well as in experimental physics. The simplest variant is a one-dimensional (Id) chain 
of coupled spin-1/2 particles, described by the Hamiltonian 

L-l 

^ = E^J^^i + ^J4m + ^K+1' (1) 
i=i 

in terms of Pauli matrices o'J''^'^ at lattice site j. Using Jordan- Wigner transformation [1] 
it can be mapped to a system of spinless fermions whose Hamiltonian has a kinetic 
(hopping) term and a density-density interaction term between fermions at neighboring 
sites. It therefore represents one of the simplest systems of strongly interacting fermions. 
The model ([1]) is exactly solvable by the Bethe ansatz [2]. Despite its solvability, 
time dependent properties, like a time-dependent current autocorrelation function that 
is via a linear response theory directly related to the transport coefficient, are at 
present beyond capabilities of exact methods. Isotropic one-dimensional Heisenberg 
model (II]) is experimentally realized in so-called spin-chain materials [3], for instance 
in many cuprates. As of yet unexplained in these materials is a very high thermal 
conductivity along the axis of Heisenberg chains believed to be due to contribution 
from Heisenberg chains and strongly influenced by impurities. 

In the present work we shall study magnetization transport in the isotropic 
Heisenberg chain; for an overview of references on a more general anisotropic Heisenberg 
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model see introductions in Refs. |5, 6J. Isotropic Heisenberg model has been studied 
in the past, however, no definite conclusion about magnetization transport has been 
reached so far. Most studies focused on the Drude weight, whose non-zero value indicates 
a non-diffusive transport, usually just called ballistic. Using the Bethe ansatz a non- 
zero Drude weight at all temperatures is advocated in [7j, with ~ 1/T behavior at high 
temperatures, while zero Drude weight at all temperatures is predicted in [8]. Exact 
diagonalizations [9] , as well as conformal field theory [10] and quantum Monte Carlo [11] , 
also result in a finite Drude weight at high temperatures, see also [121 1131 [E] • Using 
current autocorrelation function obtained by exact diagonalization is non-conclusive [T5] 
because long time scales exist. Extrapolation to the thermodynamical limit is very 
difficult with these results because exact diagonalization is limited to small systems, 
while quantum Monte Carlo, Bethe ansatz and conformal approaches have their own 
problems [6j. Low-energy bosonization calculation, together with the analysis of current 
autocorrelation function from density matrix renormalization group method, on the 
other hand indicates [6] a presence of diffusive contribution, also seen in quantum 
Monte Carlo calculation [16]. All these results, some showing zero Drude weight, 
others non-zero, or even the presence of a diffusive component, call for a more precise 
characterization of transport. Quantification just in terms of zero or non-zero Drude 
weight namely fails to distinguish between different variants of non-diffusive behavior. 

A more detailed classification can be done for instance by studying how the current 
j scales with system length L if one fixes the difference of potentials at boundaries. 
Two extreme examples are the scaling j ~ 1/L if a system obeys Fourier law, and 
j ~ in case of ballistic transport. However, as is well known from studies of classical 
transport [17], in many systems the scaling exponent is a real number, j ~ with 
< a < oo. Transport is called anomalous if the scaling exponent differs from 1, 
a ^ 1. The exponent a is via the Fourier-like law that relates a transported quantity z 
and its current j, j = —D Vz, directly related to the diffusion constant D, resulting in 
the scaling D ~ L^^". Under certain assumptions a heuristic argument with classical 
non-interacting particles leads to a connection between a and the exponent /3 of the 
spreading of disturbances as quantified by the variance a^, cr^ ~ t^. The relation is 
/3 = \lS\- The regime of 1 < /3 < 2, corresponding to < a < 1, is called super- 
diffusion, while that of < /3 < 1, corresponding to a > 1, is called sub-diffusion. Note 
that sub-ballistic transport with /3 < 2 means that the Drude weight is zero. 

Recently, numerical simulations have shown [19] that a = 0.5 in the isotropic 
Heisenberg model at an infinite temperature in the linear response regime. In the 
present work we shall extend on these results by calculating diffusion constant also 
at finite temperatures, showing that the scaling stays the same at high temperatures 
(higher than the exchange interaction). We shall also provide some other properties of 
isotropic Heisenberg model, like the relaxation rate to a nonequilibrium steady state. 
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2. The Method 

In order to be able to study nonequilibrium stationary states (NESS) we couple the 
system to reservoirs at left and right chain ends. The two reservoirs are kept at different 
potentials inducing a nonzero magnetization current through the chain. We describe 
reservoirs in an effective way using the Lindblad equation [20] for the density matrix p 
of the system, 

dp/dt = i[p,H]+C^'%p) = C{p), (2) 

where the dissipative linear operator C^^^ describing bath is expressed in terms of 
Lindblad operators L^, C'^'^'ip) = Y,k {[^kP,Ll] + [L^, pll]^. 

We use two kinds of reservoirs. To obtain NESS at infinite temperature, i.e., 
zero energy density, we use the so-called one-spin bath which is realized by two 
Lindblad operators at each end, L\ = a/1 — paf, L\ = ^yl + pa^ at the left end 
and Lf = ^/TT7^cr+, Lf = ^/T^a- at the right end, = (a^ ± icry)/2, with 0"^'^'^ 
being Pauli matrices. 

To simulate NESS at a finite temperature we use the so-called two-spin bath, in 
which one has 16 Lindblad operators acting at two boundary spins at each end. The 
form of these 16 operators is complicated and we do not state it explicitly. They are 
obtained by demanding that the stationary equation on two boundary spins, described 
by p^'^\ C'^^^{p^'^^) = 0, has for a solution a grandcanonical state p^'^\ Targeted 
grandcanonical state obtained by calculating it from a small chain of 8 spins, 

tr3^..._8(exp (— if/TL,R + /iL,RM)) (tracing is performed over 6 spins in a chain 
of length 8), where M = J2j=i'^j ^ total magnetization and Tl^r and /iL,R the 
imposed temperature and potential at the left/right end of the chain. Details of the 
implementation can be found in pi]. Note that p*^^^ is used only to generate appropriate 
Lindblad operators that will simulate finite temperature. Once a NESS is found for 
a large chain (of upto 256 spins) a real system's temperature will be determined by 
calculating expectation values of observables in the NESS. Our results do not depend 
on details of Lindblad operators used in the simulation (therefore also not on details of 
p^^-* used in deriving them), their only goal is to impose a finite energy density in the 
NESS. The driving parameter pl,r (or p) will always be small as we are interested in 
the linear response regime. For stationary properties under maximally strong one-spin 
driving see [22] . 

For both kinds of reservoirs the NESS, simply denoted by p in the rest of the 
paper, is found by evolving an arbitrary initial state p(0) with the Lindblad equation 
for sufficiently long time, until a nonequilibrium stationary state p = limi_j,oo pit) is 
reached. To calculate time evolution of p(t) we use a time- dependent density matrix 
renormalization group (tDMRG) method, as described in refs. [2T] . 

In section 13.31 we will be interested also in spectral properties of the Liouville 
superoperator £ ([2]). Because a detail knowledge of the eigenvalues of C can not be 
obtained by a simple implementation of tDMRG we use, we have instead used an exact 
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diagonalization on somewhat smaller systems. 
3. Results 

3.1. Determining temperature 

In all our simulations, using one-spin or two-spin baths, we use a weak driving 
/iL = — A* = —0.02 and /xr = +/i = 0.02. Symmetric driving with respect to left/right 
end imposes a NESS state with the average magnetization being zero. In a two-spin 
bath, with which we can simulate finite temperature states, the imposed temperature is 
the same at the left and the right end, Tl = Tr = Timp.- A consequence of this is that 
the energy density in the NESS, 

^. = (^K+i + ^^iVi + ^H+i). (3) 

is independent of the position index j and the energy current in the NESS is therefore 
zero. denotes the expectation value in the NESS, (A) = trpA. Because driving 
is weak the NESS is locally close to equilibrium. We can therefore determine [23] 
the temperature of NESS, called a "measured" temperature Tmeas., by equating the 
expectation value of the energy density to the one expected in a canonical state, hj = h^, 

^ tr (pTff) ^ exp (-g/r^eas.) . . 

^ L-l ' trexp(-i7/T„)' 
and solving for Tmeas.- The NESS can be approximated by a canonical state because the 
average potential /imeas. is zero. We have checked that the NESS state is indeed close to 
the canonical one by verifying that expectation values of nearest-neighbor observables 
in the bulk of the NESS agree with the canonical ones within (9(/iL R). For instance, for 
the NESS with Tmeas. ~ 4.8 we have one-point expectation values {cr^'^) < 10~^ (they 
are non-zero due to tDMRG truncation error), (cr|) ~ (9(/iL,R), two-point expectation 
values (cT^cT^+i) ~ ~ Wl^^l+i) ~ ^t/3, {jk) ~ C'(/iL,R), while other two-point 

nearest-neighbor expectation values are smaller than 10~^. In the thermal canonical 
state all one-point expectation values are zero, while two-point nearest neighbor agree 
with the ones in the NESS. Three-point expectation values in the NESS are all smaller 
than 10~^, except {o'k'^k+i'^k+2) well as that of all permutations of 

the three Pauli matrices occurring in these two operators) which are of order 10~^. In the 
canonical state all three-point nearest-neighbor expectation values are zero. Some four- 
point nearest neighbor expectation values in the NESS are non-zero and of size ~ 0.05, 
for instance of operators like crM+icrl+i^^l+s or ci^cr^+i^+sO^fc+s (and 

permutations of these four operators). The corresponding canonical expectation values 
are also non-zero and within 0(/iL.R) of the NESS expectation values. All canonical 
expectation values mentioned have been calculated at high temperatures by an exact 
diagonalization of small systems while the imaginary-time tDMRG method has been 
used at smaller temperatures. 

In Fig. [T]we show the dependence of the canonical energy density on temperature. 
It turns out that the boundary effects with our two-spin bath (see later) get increasingly 
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Figure 1. Average energy density in a canonical state (jlj of the isotropic one- 
dimensional Heisenberg model ([1]). Horizontal lines show energy densities hj in the 
NESS states used in the paper. 



stronger with lowering the imposed temperature. Because the operator entanglement of 
the NESS p also increases, simulations get increasingly more difficult. We are therefore 
not able to reach very low temperatures. The temperatures we used are listed in 
Table [TJ Note that the minimal temperature achieved, Tmeas. = 2.6, would in the spin 



T 

imp. 


hj 


T 

meas. 


OO 


0.0004 


OO 


50 


-0.026 


120 


10 


-0.141 


22 


5 


-0.30 


10.6 


2 


-0.69 


4.8 


0.2 


-1.17 


2.6 



Table 1. Data for NESS states used in the paper. For Ti^p = oo we use a one-spin 
bath, for others a two-spin bath. The measured temperature is determined by equating 
energy density hj in the NESS to the canonical expectation value /it- A nonzero hj 
for Tinip = OO is due to truncation errors of the tDMRG method. 

notation where H = J2j=i + + wi^h s^'^'^ = correspond 

to Tmeas. = 0.65. The exchange interaction in SrCu02 is approximately J/k^ ~ 2000 
K. Temperatures in the experiments which are of the order ~ 100 K, therefore 
correspond to dimensionless temperature T^eas. ~ 0.2 in our Pauli notation. Such low 
temperatures are unfortunately not reachable with our reservoirs [23] . 
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3.2. Magnetization profiles and the current 

The main quantity we consider is the scahng of the expectation value of the 
magnetization current in the NESS with the system size L. The magnetization current 
operator is, 

Jfc = 2Ka^+i -a^+i), (5) 

and we denote its expectation value [21] in the NESS, which is independent of the site 
index k, simply by j = (jk)- In Fig. [2] we show results at various temperatures, all for 
the same driving /xl.r = ±0.02. In addition to NESS results obtained by tDMRG we 
also show the ones obtained by numerically exactly solving [25] the master equation 
([2]). The reason that in general the current j decreases with decreasing temperature is 

0.1 



0.01 



0.001 

10 100 200 

L 

Figure 2. Scaling of current j with L for different temperatures Tmoas., all at the 
same driving /il.r- Circles are using numerically exact NESS state while squares are 
obtained using tDMRG. Data at Tmoas. = oo is the same as in Ref. [Ej, apart from 
new data point for L = 256. Two straight lines overlapping with Tmcas. — 4.8 and 
Tineas. — oo data are ^ l/L'-'-^. 

also a consequence of increasing boundary resistances due to two-spin reservoirs used. 
There is a magnetization jump at the boundary so that the first and the last spin 
have magnetization smaller than the imposed fi = 0.02. This can be seen in Fig. [3] 
for data at finite temperatures. The magnetization profiles at all temperatures from 
Table [1] (except the ones at Tmeas = 2.6) can be well described by the scaling function 
(o"^) = fc^ arcsin(x), where x = 2^-^^ — 1 is a scaled position and A; is a temperature 
dependent prefactor, effectively taking into account for boundary jumps. For instance, 
at Tmeas. = 10.6 (left frame in Fig. [3]) it is /c ~ 0.65. Because of the boundary jumps, to 
properly account for the scaling of j with L, ie., to asses the validity of the Fourier law, 

j = -DV,(0, (6) 
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Figure 3. (Color online) Left: Magnetization profiles at Tmcas. = oo and at 
Tmcas. = 10.6. At lower temperature the magnetization exhibits jumps at the boundary. 
At all temperatures the profile is well described by the scaling function arcsin(a;). 
Right: At lower Tmcas. = 4.8 the convergence of profiles to arcsin (x) seems to happen 
at larger sizes L than at higher temperatures. 



we have to scale the current with the actual magnetization difference given by 
(al) — (a^) = Az. This is shown in Fig. HI From the figure we can estimate that 



0.01 



"^meas. 


=2.6,4.8,10.6,22,120,00 
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1.55/L°-^ 
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Figure 4. Current divided by Az = {af — a\) for different L and temperatures. For 
sufficiently large L all seem to converge to the same line ~ 1.55/L°'^. 



the current is 
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independent of the temperature (energy density), at least for Tj^eas. ^5. In this range 
of temperatures the diffusion constant D ([6]) therefore scales as 

L)^1.55L°-l (8) 

At lower temperatures it is difficult to assess if the scaling is still the same. Data in 
Fig. Hlfor Tmeas. = 2.6 and 4.8 show larger current than predicted by Eq.©, however, 
one explanation could be that the length at which the asymptotic behavior ([71) begins 
gets larger than the sizes studied. In the right frame of Fig. [3] we can for instance see 
that the asymptotic arcsinx magnetization profile is at Tmeas. = 4.8 not yet reached 
for L = 32, whereas at higher temperatures this happens already for smaller L's (left 
frame) . 

Scaling of the current j ~ 1/y/L or, equivalently, of diffusion constant D ~ \/L, 
can be used to show in a non-rigorous way the spatial dependence of the magnetization. 
We shall use the Fourier law (EI) with a space-dependent diffusion constant D{r) at 
site r. Because equation ([8]) tells us that the diffusion constant gets larger for larger 
chains, close to boundaries, local diffusion constant should become smaller. Assuming 
a square-root scaling we must have D{r) oc \/r[L — r)/ L. This gives a differential 
equation 

const. r(L — r)dz , . 

where z = (cr^). Integrating the above equation with appropriate boundary conditions 
one immediately gets the profile z ~ arcsin {2r/L — 1). 

3.3. Spectral properties of C 

Spectral properties of a Liouvillean superoperator £, i.e., the linear operator 
representing the right-hand-side of the master equation ([2]), are important for several 
reasons. For instance, they determine the relaxation rate to NESS as well as deviations 
from NESS expectation values at finite times. The superoperator C is non-Hermitean 
and therefore has a spectrum of eigenvalues A^, k = 0, ... ,4^ — 1, lying in a complex 
plane. We shall order eigenvalues Xk in a descending order according to their real part, 
starting with the largest Aq = 0. The right eigenvector |xq ), corresponding to Aq, is the 
sought-for NESS state p, symbolically \ x^) =\ p). As a consequence of trace preservation 
of C the left eigenvector (xq | corresponding to Aq is on the other hand proportional to 
the identity operator ~ 1, irrespective of the system, (xqI = In our system Aq is 
always nondegenerate. For simplicity we shall in this subsection discuss properties of 
the isotropic Heisenberg model with a one-spin bath, that is at an infinite temperature. 
The driving potential is weak, = 0.02, however, the values of the eigenvalues are in 
the linear response regime largely independent of p. 

Besides eigenvalues and eigenvectors that are closest to Aq (in real part) are 
also of interest. Because the dynamics governed by the Lindblad equation is contractive 
all real parts of eigenvalues are non-positive, Re(Afc) < 0. Linear operator C is in 
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general non-diagonalizable with the Jordan canonical form, see for instance Ref. [26] for 
a discussion of spectral decomposition in such case for quadratic fermionic systems. For 
isotropic Heisenberg model ((21) we have found by numerical computation that for few 
eigenvalues with the largest real parts there are always as many linearly independent 
eigenvectors as is the multiplicity of the corresponding eigenvalue and therefore the 
Jordan form for these eigenvalues is trivial of dimension 1. We can therefore write 

C = \o\x^) {x^ I + Ai I xf) {x\ I + Aa I x^) (x^ I + ■ ■ • , (10) 

where left and right eigenvectors are mutually orthogonal, = Sjk, with the 

standard Hilbert-Schmidt inner product {A\B) = tr(A^S), and we normalize left 
eigenvectors {x^ |. Then the value of the real part of Ai, also called the gap of Liouvillean, 
determines the convergence rate with which the NESS is reached from p(0), while the 
corresponding eigenvector gives the deviation of p{t) from the NESS p. In addition, 
the scaling of the gap with the system size L can be used to locate nonequilibrium 
phase transitions. For studies of this phenomenon in quantum system see Ref. [27], for 
classical systems see e.g. [28] . 

We have determined the lowest 4 eigenvalues Ao,i,2,3 and their corresponding left 
and right eigenvectors using numerically exact diagonalization on small systems of size 
L < 12. In addition, we determined the relaxation rate r of our tDMRG solution p{t) 
by fitting the convergence of magnetization at the middle of chain to its asymptotic 
value, tr {p{t)a^/2) ~ ~ ^^P (^^^ same r is also obtained by looking at the 

convergence of magnetization current). In our isotropic Heisenberg chain with one-spin 




Figure 5. Gap of the superoperator £, eq. ©, for a one-spin bath (i.e., infinite 
temperature) from an exact calculation (squares and circles) and relaxation rate of 
tDMRG simulations r (crosses). The largest nontrivial eigenvalue of £ scales as 

Ai lOO/L^" (Ml line), while the 2nd one goes as A2 (dashed 

line). The scaling of tDMRG convergence rate r with L is the same as that of A2. 
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bath Ai and A2 always have zero imaginary part, Im(Ai^2) = 0. Data in Fig. Elshows that 
the gap of C decreases as Ai ~ The same scahng with L is obtained also for the 

XX model [29]. The eigenvalue Ai is 2x degenerate for L < 6, while it is nondegenerate 
for L > Q. On the other hand, A2 becomes 2x degenerate for L > Q. The scaling of A2 
for large L is A2 ~ and therefore decays with L in a much slower way than Ai. 

What is interesting is that the convergence rate of tDMRG simulation is not given by Ai, 
but rather follows the scaling of A2. This is so because expectation values of current and 
magnetization, which are relevant observables for our discussion of transport, are very 
small in the eigenvector corresponding to Ai. In fact, for L < 6 they are identically zero, 
while for larger L their values are shown in Fig. El One can see that the contribution 




Figure 6. The current in the NESS state, j (squares; the same data as in Fig. [2]), 
and the current in the 2nd eigenvector, corresponding to Ai, at the 1st site (where it is 
the largest), = {ji\xf^). All using exact diagonalization and a one-spin bath (i.e., 
T = 00). 

from I xf ) to the magnetization current (and magnetization as well) scales as ~ 1/L^ and 
is indeed negligible in the thermodynamic limit |30]. For instance, relative contribution 
at L = 128 would be (0.18/L2)/(0.062/L°-5) ^ 0.002, which is below the precision of 
our tDMRG simulations. Namely, we estimate that the truncation errors in tDMRG 
simulations result in relative error in the current below 1% at Tmeas. = 00 and below 5% 
at Tjneas. = 4.8, both for the largest sizes shown. 



4. Conclusion 



Using extensive numerical calculations of nonequilibrium steady states close to 
equilibrium in chains of upto 256 spins we have shown that the diffusion constant 
of magnetization in the isotropic Heisenberg model scales with the system length as 
D ~ L"'^ for temperatures larger than the value of the exchange interaction. In this 
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temperature regime the anomalous diffusion exponent of 0.5 seems largely independent 
of the temperature. The spectral properties of the Liouville superoperator have also 
been explored, showing that that the gap scales as ~ with the system length L, 

while the tDMRG expectation values of magnetization and current converge on a faster 
time increasing as ~ L^'^^. 
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